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We study formation of black holes in the Friedmann universe. We present a formulation of the 
Einstein equations under the constant mean curvature time-slicing condition. Our formalism not 
only gives us the analytic solution of the perturbation equations for non-linear density and metric 
fluctuations on superhorizon scales, but also allows us to carry out a numerical relativity simulation 
for black hole formation after the scale of the density fluctuations is well within the Hubble horizon 
scale. We perform a numerical simulation of spherically symmetric black hole formation in the 
radiation-dominated, spatially flat background universe for a realistic initial condition supplied 
from the analytic solution. It is found that the initial metric perturbation has to be non-linear (the 
maximum value of 3D conformal factor ifio at t = should be larger than ~ 1.4) for a black hole 
to be formed, but the threshold amplitude for black hole formation and the final black hole mass 
considerably depend on the initial density (or metric) profile of the perturbation: The threshold 
value of t[>o at t = for formation of a black hole is smaller for a high density peak surrounded by a 
low density region than for that surrounded by the average density region of the flat universe. This 
suggests that it is necessary to take into account the spatial correlation of density fluctuations in 
the study of primordial black hole formation. 

PACS: 04.25.Dm, 95.35.+d, 97.60.Lf 



I. INTRODUCTION 

The formation of black holes in the early universe and 
its cosmological implications have been discussed in a 
variety of contexts for decades |Q. However, it has long 
been thought that it would be practically impossible to 
prove or disprove the existence of these primordial black 
holes. Recent discoveries of microlensing events by MA- 
CHOs of mass ~ 0.5M Q in the halo of our galaxy [g] have 
dramatically changed this situation. By various other 
means of observation, it may be possible to deny all the 
other possibilities and hence to identify MACHOs with 
black holes. Then these MACHO black holes must be 
primordial since it is impossible to form a black hole of 
mass smaller than ~ Mq as a result of stellar evolution 
Furthermore, it has been recently suggested that if 
MACHOs are in fact primordial black holes, the number 
of binaries that are just coalescing today may be large 
enough to be directly detected by the oncoming gravita- 
tional wave observatories such as LIGO, VIRGO, GEO 
and TAMA within a few years |Q. Consequently, it has 
become an urgent issue to quantify how and when these 
black holes could be formed in a precise manner. 

Among other possibilities, primordial black holes are 
most conceivably formed from large curvature perturba- 
tions generated during an inflationary stage of the very 
early universe [js],^)- The curvature perturbations gener- 
ated in the inflationary universe are dominated by the 



so-called growing adiabatic mode of density perturba- 
tions. In the linear theory, the evolution of these per- 
turbations is well studied and their temporal behavior 
is known throughout the whole stage from the epoch 
when their wavelengths are much larger than the Hubble 
horizon scale until their evolution becomes non-linear on 
scales much smaller than the Hubble horizon. 

However, the amplitude must be already large enough 
(of order unity) to form black holes when the character- 
istic wavelengths of the perturbations were on superhori- 
zon scales. Furthermore, the formation of a black hole is 
itself a fully general relativistic phenomenon. The evolu- 
tion of non-linear density perturbations on superhorizon 
scales was investigated by several authors and the thresh- 
old amplitude of curvature perturbations on superhorizon 
scales for forming black holes was estimated . These 
previous estimates of the threshold amplitude were based 
on approximate analytical treatments and/or on rather 
naive numerical simulations, hence are admittedly crude. 
Furthermore, there is a crucial reason that requires us an 
accurate estimate as follows: According to the inflation- 
ary scenario, the probability distribution of the curva- 
ture perturbations is essentially Gaussian and primordial 
black holes are produced from the high amplitude tail of 
the distribution. Therefore, a small error in the estimate 
of the threshold amplitude will result in a large error in 
that of the number of produced black holes. Thus thresh- 
old amplitude must be estimated accurately. 
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As a first step to accomplish this purpose, in this pa- 
per, we present a new formalism by which it is possible to 
follow the formation of a primordial black hole through- 
out the whole stage starting from the very early uni- 
verse when the perturbation is well outside the Hubble 
horizon to the final stage when a black hole is formed. 
More specifically, our formalism not only gives us the 
analytic solution of non-linear curvature perturbations 
on superhorizon scales but also allows us to perform a 
numerical simulation of the black hole formation with 
the initial data given by the analytic solution, with no 
need of changing the gauge conditions or of numerical 
matching. In addition, it may be worthwhile to mention 
that the constant mean curvature time-slicing employed 
here is equivalent to the so-called constant Hubble slic- 
ing in cosmological perturbation theory |)| . And, it has 
been pointed out that the constant Hubble slicing is most 
appropriate for evaluating non-linear curvature fluctua- 
tions generated during inflation jl0|]. Hence the initial 
curvature perturbation spectrum evaluated in models of 
inflation can be directly used for the initial data of our 
problem. 

Then using our formalism, we carry out a spheri- 
cally symmetric simulation of black hole formation in 
the radiation-dominated Friedmann universe. We con- 
sider the initial data with two parameters; one describes 
the amplitude and the other the radial profile. We find 
that both the threshold amplitude for black hole forma- 
tion and the final black hole mass depend appreciably 
on the initial profile of the perturbation. We also con- 
sider another possible criterion for black hole formation 
by defining a compaction function of the perturbation 
profile. Although this function can be defined only for 
a spherically symmetric configuration, we find the maxi- 
mum value of this function gives us a better criterion for 
the formation of black holes. 

While this paper was in preparation, a couple of pa- 
pers on the primordial black hole formation appeared 
on the astro-ph fu]| . It seems that their formalism is 
powerful for studying the formation of a black hole in 
a spherical symmetric spacetime. However, it does not 
seem convenient to give a realistic initial condition which 
should be supplied just after inflation. Actually, they 
give initial conditions at the epoch when the scale of 
the density fluctuation is as small as the Hubble horizon 
scale. Since the density fluctuation is already nonlinear 
at that epoch, it is impossible to control the initial date 
so that it reduces to the growing adiabatic mode when 
the evolution is traced back in time to the very early 
universe. In other words, their initial conditions are in- 
evitably contaminated by unrealistic decaying mode per- 
turbations which badly diverge as t — > 0. As a result, 
though the criterion for black hole formation they find is 
new and interesting, it cannot be directly related to the 
initial condition at the end of inflation. So that it is not 
convenient for a practical study of primordial black hole 
formation. Furthermore, application of their formalism 
is restricted to the spherical symmetric case (i.e., very 



special case). In contrast, in ours, it is easy to relate a 
criterion of black hole formation to an initial condition 
just after inflation, and also it can be applied to general 
3D cases. The only restriction of our formalism is that 
the spacetime be asymptotically spatially flat Friedmann. 

The paper is organized as follows. In Sec. II, we 
present the Einstein and hydrodynamic equations in the 
Friedmann universe using the 3+1 formalism, which have 
appropriate forms for numerical relativity simulations. 
We then introduce the constant mean curvature time- 
slicing in which the equations for geometric variables 
have similar forms to those in the maximal slice con- 
dition in the asymptotically flat spacetime. In Sec. Ill, 
assuming that the length scale of a density fluctuation 
is always much longer than the Hubble horizon scale, we 
take the long wavelength limit of the equations derived 
in Sec. II, and then find the analytic solution for the per- 
turbation equations. In Sec. IV, we perform numerical 
simulations of black hole formation in a spherically sym- 
metric, radiation-dominated universe using initial condi- 
tions given by the analytic solution in Sec. III. Sec. V is 
devoted to summary. Throughout this paper, we use the 
units c = 1 = G. 



II. FORMULATION 

We present the Einstein and hydrodynamic equations 
in the Friedmann universe using the 3+1 formalism in 
general relativity. We write the line element as 

ds 2 = g^dx^dx" 

= {-a 2 + /3 k (3 k )dt 2 + 2(3idx l dt + -..d.r'd.r 1 . (2.1) 

where g^, a, (3 l ((3i — 7ij/? J ), and 7ij are the 4D metric, 
lapse function, shift vector, and 3D spatial metric, re- 
spectively. Since we consider an asymptotically spatially 
flat Friedmann universe, we rewrite 7^ as 



7y = ip A a(t) 2 7ij, 



(2.2) 



and impose the condition det(7y ) = Aet{rjij) = rj, where 
rjij is a flat spatial metric. Here, a(t) is defined to be the 
scale factor in the homogeneous universe, i.e., the scale 
factor in the asymptotic region, and we determine it from 
the well known equations for the scale factor as 



4-7T 



Stt 



p (t) + 3P (t) 



= -5-0 Po(t), 



(2.3) 
(2.4) 



where po and Pq are the density and pressure for the ho- 
mogeneous universe, and d = dta. As we find below, ao, 
po and Po are automatically determined when an equa- 
tion of state is provided. 



We also rewrite the extrinsic curvature as 



Kij = 1/1*0? A. 



(2.5) 
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where K = K k and hence Aij is defined to be traceless. 
The indices of and A y are to be raised or lowered in 
terms of 7^ and 7^. In numerical computation, we will 
solve 7^, A^, if] and K instead of 7^ and K^j. Here- 
after, we use Di and Di as the covariant derivatives with 
respect to 7^ and 7^ , respectively. 

As a source of the energy momentum tensor, we con- 
sider a perfect fluid for which the energy momentum ten- 
sor is written as 



(p + P)u^u v + Pg^, 



where p, P and are the energy density, pressure and 
four velocity, respectively. Hereafter, we assume a rela- 
tivistic polytropic equation of state, P = (T — l)p, where 
T is a polytropic constant and we assume T > 1 in the 
following discussion. For the relativistic polytropic equa- 
tion of state, we get the following relations from Eqs. 

& 

a = a f t 2 ' 3T and p = ^5^, (2.7) 

where aj is a constant. 

The hydrodynamic equations are written in the form 

d t (w^a 3 p^ r ) + -l-d k (r, 1 / 2 w^ (i a 3 p 1/r v k ) = 0, (2.8) 
d t {wip 6 a s (p + P) Uj } + -La^Va^Vfc + P)v k Uj} 
= -a^ 6 a 3 d 3 P + w^a 3 (p + P){-au°d j a 

Wi^-I^^}, (2.9) 

where w = au°, and 



~M u i 



(2.10) 



Evolution equations for geometric variables are written 
as follows [0 : 



(d t - (3 k d k )% = -2ai ij 

2 



(d t - p k d k )A iJ 



(2.11) 



- ( DiDjOt - ^-D k D k a 



+a(KAij - 2A lk A j 



+f3 k ^A kj + k tj A ki - -^D^An 



_8_7ra . 



(-•' :-! 

4> 



7ij g k 



(d t - f3 k d k )4> + ^ = i{~ aK + (O)D0 k } 
(dt - p k d k )K = a(A l0 A^ + ^K 2 ) 

-D k D k a + ATta(E + S k k ), 



(2.12) 
(2.13) 



where Rij is the Ricci tensor with respect to 7^ , (q) D k 
is the covariant derivative with respect to rjij , and 



(p + P)Ui1lj + Pji L 



(2.15) 



To clarify the meaning of equation (2.12), we rewrite 



Rij as 



R t j - Rij +Rij, 



(2.16) 



(2.6) where Rij is the Ricci tensor with respect to 7y and 



6 - 2 ~ u 

+—DiipDjij - —%D k ^D k 4>. 



Rij is written as 



(2.17) 



Ri 



■(o)A7y + (o)Dj(o)D k j ki + (o) A(o)-D K 7fei 



(2.18) 



where (o)A is the Laplacian with respect to rjij, f 



k! 



~ hi 

\ ((o) A7iJ + (0)-Di7« ~ (o) A7y) ■ (2.19) 



Note that we use a relation 7^ (o yDfc7 ij = rf"* (QjDkVij = 
to derive Eq. fl2.18| ). From Eq. (|t|), it is found that un- 
der an appropriate gauge condition such as a transverse- 
t racel ess (TT) gauge, ^D k ^f kl = 0, Eqs. ( [2.11 ) and 
(|2.12 ) are found to constitute a wave equation for tensor 
lij- 

Hamiltonian and momentum constraint equations are 



R k k - AijA ij + -K 2 = 16irE, 



DiA l j - -DjK = 8nJ 3 , 
o 



where 



E = (p + P)w 2 - P, 
Ji = (p + P)wui, 

We may write the constraint equations as 
Aip = ^-tp - 2n^ 5 a 2 E 



(2.20) 
(2.21) 



(2.22) 
(2.23) 



-K' 



&(iP 6 A tJ ) - -iP 6 D t K = 8ttJ^ 6 , 



(2.24) 
(2.25) 



(2.14) where A is the Laplacian with respect to 7^. 
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In this paper, we choose a constant mean curvature 
slice as we chose in a previous paper Ea] 



K = K(t) = - 



3« 



(2.26) 



This choice can most effectively factor out overall factors 
of the expansion of the background universe from the 
dynamical variables. In this case, we obtain the equation 
for a as 

Aa = a[47r{2(p + P){w 2 - 1) + p - Po + 3(P - P Q )} 



+ A ij A ij ] + 12Tr(p + P Q )(a-l), 



(2.27) 



where A is the Lap lacian with respect to 7, j . We also 
note that Eqs. ( 2.1 3| ) and ( [2.24| ) are, respectively, rewrit- 
ten as 



-±^(^(3%, (2.28) 



(d t - (3 k d k )^ = -ip{a - 1) 
A0> = ^0 " 27n/> 5 a 2 [( j0 + P){w 2 - 1) + p - po] 

O l J 



(2.29) 



Thus, in the constant mean curvature slice condition, the 
equations for a and ip are similar to those for the max- 
imal slice condition, K = 0, in the asymptotically flat 
spacetime. Hence, we expect that it has a singularity 
avoidance property for the case of black hole formation. 



III. LONG WAVELENGTH LIMIT 

In this section, we consider the so-called long wave- 
length approximation assuming that the characteristic 
length scale L of a density fluctuation is always much 
larger than the Hubble horizon scale, L ^ t ^ a/a. 
First, we introduce a small parameter e, and assume that 
S = p/ pa — 1 is of 0(e 2 ) and its characteristic length scale 
is of 0(l/e). The latter assumption is equivalent to as- 
suming that the magnitude of spatial gradients of the 
quantities is given by diip — ip x O(e), dia = a x O(e), 
did = S x 0(e) and so on. Then, it is found from the 
equations presented in Sec. II that the following relation 
should hold: 

^-l = O(e ), 
u* = 0{e l ), 

Aij, hij(= 7y - %■), x(= " ~ 1), 8 = 0(e 2 ), 

Ui , v l +fi l = 0(e 3 ). (3.1) 

Here we have assumed for simplicity that the amplitude 
of primordial gravitational wave perturbations is negli- 
gible. We note that because the gravitational wave per- 
turbations do not decay with time evolution, the follow- 
ing solution can change considerably if their amplitude is 
not small initially. We have also assumed that the vector 



part is absent for any quantity; in other words, we do not 
consider vorticity. 

It is worth mentioning that these assumptions are nat- 
urally realized in most of successful inflation models. 
In the inflationary universe scenario, only the so-called 
growing mode perturbations of scalar and tensor types 
survive and amplitude of the tensor perturbation is gen- 
erally very small. We note that -0 — 1 may be larger than 
unity, i.e., we have not used any approximation for treat- 
ing ip. On the other hand, other quantities should be 
small enough and are regarded as small perturbations. If 

— (i.e., if the linear theory applies), the above 

corresponds to assuming the existence of only the adia- 
batic growing mode. 

Because we have not yet imposed any condition on (3 k , 
we cannot specify the order of magnitude of (3 k and v k . 
For example, in the case of the minimum distortion gauge 
Q , they are of 0(e). On the other hand, in the k = 
gauge, v k = 0(e 3 ). In the following discussion, we do 
not have to specify the gauge condition for (3 k in order 
to obtain solutions except for (3 k and v k . 

Substituting the lowe st o rder terms in O(e) of each 
variables shown in Eqs . (3.1), we have the following equa- 
tions: 



h+^ + V k v k = 0(e% 
1 ip 

d t {p a 5 V k (v k + (3 k )} + ^ Po a 5 V k (v k + (3 k ) 



(3.2) 



+Poa V fe 



(3 k )} 



= - P0 « 3 (V 2 x+^V 2 <5j+O( e 4 

6 -i-Z-X = V k (3 k + 0(e% 
ip a 

(o) A0 = -2™VVo£ + O(e 4 ), 

V 2 x = 47rp a 2 {(3r - 2)6 + 3T X } + 0(e 6 ), 

2 
3' 



d t hij = -21^ + Sik/Pj + 5 3k f3 k tl - -SijP k k 



(3.3) 

(3.4) 

(3.5) 
(3.6) 

(3.7) 



a 



1 



ip 4 a 2 



((o)A(o)-Dj^ - -y (o)Ai/>) 



+ ~Jp. {(o)Diil}(o)Djil> - —( )D k ip( )D k ip J 



iP 2 V"> "-w-^ 3 
+ 0(e 4 ), (3.8) 



where 



v 2 = 



(3.9) 



4 



The first two equations, (3.2) and (p.3|), a re d eri ved from 



hydrodynamic equations and other five, (3.4) — (3 
derived from e qua tions f or g eometric variables 



From Eqs. (3.2) and (3.4), we find that the following 
relations have to hold; 

^ = 0(e 2 ), (3.10) 

h + 3- X = -W k (v k + (3 k ) = 0(e 4 ). (3.11) 
1 a 



Also, we find that the right-hand side of Eq. (3.6) has to 
be of 0(e 4 ), i.e., 



(3r - 2)6 + 3T X = 0(e 4 



(3.12) 



From these relations and a reasonable assumption that 
5^0 for t — > 0, we finally find the time dependence for 
each variable at the leading order in O(e) as 



3r 



:Xoct 2 - 4 / 3r , 



3r- r 

V fe (« fc +/5 fe )oc^ 8 / 3r , 

Wfe °ci 3 - 4 / 3r , 
ip oc t°, 

An cx t 1 -^. 



(3.13) 

(3.14) 
(3.15) 
(3.16) 
(3.17) 



Time dependence of v l , /3\ 0(e 2 ) part of tp, and hij is 
found when we give spatial gauge condition for (J 1 . For 
example, in the TT gauge and/or minimum distortion 
gauge fll4| , it is found that 



v k , [3 k cx ^- 4 / 3r , 

hij oc i 2 " 4 / 3r , 

0(e 2 ) part of tp cx t 2 - 4/3F . 



(3.18) 
(3.19) 
(3.20) 



We emphasize again that we do not restrict our atten- 
tion to the case where tp — 1 <C 1. Namely, even when 
the scalar part of the metric is non-linear, we can still 
find the analytic solution as long as the long wavelength 
approximation holds. 

The purpose of this paper is to give a framework to 
investigate the primordial black hole formation process, 
and its standard formation scenario is as follows: In 
a very early phase of the universe, just after inflation, 
the scalar-type perturbations generated from the quan- 
tum fluctuations of an inflaton field have the length scale 
much larger than the Hubble horizon scale at that time. 
Some of these perturbations may have a large metric per- 
turbation amplitude J5|-||. As long as its scale is larger 
than the Hubble horizon scale, it never collapses, but 
once the scale becomes smaller than the Hubble horizon 
scale, it collapses to form a black hole. 

The advantage of our present formalism is as follows: 
Once we give a realistic initial condition at the very early 
epoch just after inflation, evolutions of the metric and 
density fluctuations can be analytically calculated as long 
as the length scale of the fluctuation is much larger than 



the Hubble horizon scale. Then we may start a numerical 
simulation sufficiently before the scale enters the Hubble 
horizon and follow the black hole formation process with- 
out changing the gauge condition or numerical matching 
of initial data. Hence, we can consistently investigate the 
evolution of the metric and density fluctuations through- 
out the whole dynamical range starting from a very early 
epoch of the universe at which analysis can be performed 
in an analytical manner up to the formation of a black 
hole when analysis should be done in numerical relativity. 



IV. NUMERICAL STUDY IN THE SPHERICAL 
SYMMETRIC CASE 

To demonstrate the usefulness and robustness of our 
formalism, in this section, we perform numerical simula- 
tions of primordial black hole formation assuming spher- 
ical symmetry. 



A. Basic equations 



For the spherical case, the line element can be written 



as 



ds 2 = -(a 2 - il) 4 (3 2 r 2 )dt 2 + 2il) i a 2 (3rdrdt 
+iP 4 a 2 {dr 2 +r 2 dn), 



(4.1) 



where (3 denotes (3 r jr. We may say that we choose the 
minimum distortion gauge condition in this line element 
because dtfij — |L4|]. As we mentioned in Sec. II, it is 
not always necessary to take the spatial gauge condition 
used here, and we may use other gauge conditions such 
as (3 r = 0. The spatial gauge condition we choose here is 
only one example p7| . 

Since gravitational waves are not generated in the 
spherically symmetric spacetime, we need not solve the 
evolution equations for the geometrical variables if we 
solve the constraint equations and the equations of the 
gauge condition. On the other hand, if we solve the evo- 
lution equations, we can use the constraint equations in 
order to check the accuracy of numerical solutions in each 
time step. Thus, we solve the evolution equation for tp 
instead of the Hamiltonian constraint equation, but use 
the latter to check the numerical accuracy. Then, the 
equations for the geometric variables solved in numerical 
computation are as follows: 

(d t - 2(3yd y )i> = ^(a - 1) + |(3/3 + 2yd v {3), (4.2) 

C0) AC = CV> 4 a 2 [6tt(p + P)(w 2 -1) 

+2tt{ p - p Q + 6(P-P )} 



-12n(p + P ) 



21 
16' 



-A" 



+ip 5 a 2 4n{2(p + P){w 2 - 1) 
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+p - po + 3(P - P )} + 



2/9 y /3 = -aA, 
d r (r 3 iP 6 A) = 87rr 4 V> 6 (p + 



(4.3) 

(4.4) 
(4.5) 



where y = r 2 , u — u r /r and £ = ip(a — 1). We use the 
relation -A g e /2 = -I//2 = A r r = A. Eq. @) is 
solved under the boundary condition at r — > oo as 



c 



H e -r/r 



0(r~ 



(4.6) 



where C is a constant, and ro = 1 / ^/ 12ttT poa 2 . 
The Hamiltonian constraint equation 



(o)Ai/> 



-2Tnp 5 a 2 [{p + P){w 2 
3ip 5 a 2 



1) + P - H 



1G 



(4.7) 



is solved only on the initial time slice with the outer 
boundary condition ip — > 1 + C^/2r + 0(r~ 3 ) where C^, 
is a constant. 

The hydrodynamic equations are written in the form 

d t (w^ 6 a 3 p 1/r ) + r- 2 d r (r 3 wij 6 a 3 p^ r v) = 0, (4.8) 
d t {wfa 3 (p + P)u) + r- 3 d r {r 4 w^a 3 {p + P)vu) 

= -2aiP 6 a 3 d y P + wip 6 a 3 (p + P)i-2wd y a 

Aau 2 y 



(4.9) 



wher e v = v r /r. Using the relation, p 



^constant, 



Eq. (4.8) may be written as 

d t (wifj e e) + r- 2 d r (r 3 wip e ev) = 0, 



(4.10) 



where e = (p/po) 1 ■ In numerical simulation ,we choose 
D = wip 6 e and S = wip 6 (p + P)u as variables to be 
solved. Once D and S are given, w is obtained by solving 
the algebraic equation 



p 2 ^ 2 T 2 



2V 



(w 2 1) 



S V_ w 2T-2_ 



ip 4 a 2 



(4.11) 



and v is then given from 

cm 



v = -/3 + 



wip 4 a 2 



= -/3- 



a5 



Tw; 2 V' 1 Voa 2 P/^ 6 ) r ' 



(4.12) 



To examine whether a black hole is formed or not, in 
each time step we search for apparent horizon which is 
defined as the outermost trapped surface |l5[ |. In the 
spherically symmetric case, the outermost zero point, r = 
tah, of the function 



Q(r) 



2 a - + A 
a 



1 

ip 2 a 



(4.13) 



corresponds to the outermost trapped surface p3[ . 
Hence, we only need to calculate 6 in each time step 
and look for the zero point. 

Once the apparent horizon is determined, we also cal- 
culate the mass of the apparent horizon which we define 
as 



A/ah 



aip 2 



at r = r AH - 



(4.14) 



Mah is not identical with the gravitational mass of a 
black hole in general. However, if it settles down to a 
constant in the late epoch after formation of the black 
hole, it can be regarded as the gravitational mass because 
in the spherical and static space Mah agrees with the 
gravitational mass. 

In order to check the numerical accuracy of Mah, we 
also calculate the conserved Kodama mass in the spher- 
ical spacetime jl6| . The Kodama mass within a radius r 
is defined as 

M K (r) = Ait [ r' 2 dr'a 3 a^T\K^, (4.15) 
Jq 



where the components of are 
1 d 



K 



K r = 



aip 2 dr 



ip 2 r, 



l_d_ 

tip 2 dt 



V> r, 



(4.16) 
(4.17) 



and K — K v — 0. Since M K at r = tah is proved to be 
equal to A/ah |@, it can be used to check the accuracy 
of our estimation of A/ah- In our simulations, we found 
that both agree well except for the very late epoch after 
formation of a black hole at which the gradient of a near 
the apparent horizon is very steep and it is difficult to 
keep numerical accuracy well. 

Numerical simulations are performed taking 3000 inho- 
mogeneous grid points for the r-axis. The circumferential 
radius of the outer boundary is always kept to be much 
larger than the Hubble horizon scale. We also take a suf- 
ficient number of grids inside a black hole forming region. 
More concretely, we take grids as = Ar(f l — 1)/(/— 1), 
where i = 0, 1, 2, • • • 3000, / is a constant slightly larger 
than 1 , and Ar is the grid spacing at origin which is much 
smaller than the circumferential length of a formed black 
hole. 



B. Initial conditions 

We make use of the analytic solution derived in Sec. 
Ill to give a realistic initial condition for p, u and so on. 
Thus, the initial condition is given at an epoch when the 
length scale of a density fluctuation is much larger than 
the Hubble horizon scale. 

First, we assume that i5 = p/po — 1 is much smaller 
than unity and write it as 
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5 = f(r)t 2 - 4 / 3r . (4.18) 
FromEq. fljTlf) , we soon get % = a — 1 as 

3r- 2 



3r 



-f(r)t 



2-4/3r 



(4.19) 



m is derived from Eq. (4.9) in the long wavelength limit 

as 



3r + 2 



(4.20) 



Thus, if we specify the function f(r) on the ini tial time 
slice t — t , and subsequently solve Eqs. (4.3)— (4.5) and 
(4.7), we obtain the initial data for A, (3, v, ip and a. 
In this paper, we simply give 6 as 



S • V> 6 = C s 



exp 



(-^- CT ~ 3cxp (-^y 



2-4/ST 



(4.21) 



-2/3T 



where Cs, ro and o are constants, and we set a/ 
Roughly speaking, Cs and a specify the amplitude and 
shape of the density fluctuation, respectively. r Q deter- 
mines the length scale of the density fluctuation, and we 
fix it to be unity in the following. Hence, hereafter, the 
mass and length are shown in the units ro = 1. 

If we define the spectrum of density fluctuations as 



S(k) 



jo(kr)8 ■ ip r dr, 



(4.22) 



where jo(kr) = sm(kr)/kr is the spherical Bessel func- 
tion of the 0-th order, we get 



6(h)-- 



t 

ro 
-k 2 r, 



2-4/3T 



C s 



° XP ( — 4~^) ~~ 6XP (" 



(4.23) 



Thus the wavelength of the dominant spectral compo- 
nents of the density fluctuation is larger than ~ nr in 
the comoving scale (or k < 2/ro). Hence, we start all the 
simulations at an initial time t which satisfies the condi- 
tion t -C 7r(l — 2/3r)a(t)ro- In the following, we always 
set the initial time as 10 _4 ro. 

Initial conditions are numerically determined by per- 
forming iteration as follows: (a) We so lve E q. ( |4.7| ) for 
ip, for the density profile given by Eq. (4.21). (b) From 
S ■ ip 6 as well as ip obtained in (a), we determine f(r) . 
(c) we calculate u, w, A, and /3 by using Eqs. ( 4.20 ), 
w = y/\ + u 2 r 2 ip~ 4 , (4.5) and ([hj), resp ectively, and 



substitute the new w and A into Eq. (4.7). We repeat 
this procedure until sufficient convergence is achieved. 

Hereafter, we pay attention only to the r = 4/3 
case, because it is most probable that the universe was 



radiation-dominated in the early times. In this case, for 
t — ► 0, the metric behaves as 



[3 — > const <C 1, 

ip — 1 — *■ const = 0(1). 



(4.24) 
(4.25) 
(4.26) 



Apparently, the metric is non-linear for large Cs at t = 
because ip - 1 = 0(1). Note that 5 -> 0, for t -> 0. 
Hence, in the present coordinate condition, a black hole 
is formed only if the metric is non-linear even though the 
density fluctuation is very small. In other word, we may 
say that the criterion of formation of black holes depends 
only on ip at t = 0. 



C. Numerical results 

Numerical simulations were performed for various val- 
ues of Cs and a. Specifically, we surveyed a two- 
dimensional parameter space with Cs > and <j > 1. 
We note that in the case a = oo, the high density peak 
is surrounded by a flat universe, and in other cases, it 
is surrounded by a region in which the density is lower 
than that of the flat universe. 

We show the spectrum of 5(k) for the cases a = 1.5, 
2, 3, 5, 8 and oo in Fig. 1. For a — oo, the peak is 
at k = 0, while in the other cases, the peak wavenumbcr 
and the width around the peak becomes larger and wider, 
respectively, for smaller a. 

In Fig. 2, we show ipo — ip(r — 0) at t — as a function 
of Cs for a = oo, 2, 3 and 5. Here, we can easily calculate 
■00 at t — by using the analytic solution found in Sec. 
III. It is found that ipo at t — monotonically increases 
with increasing Cs, irrespective of a. 

In Figs. 3 and 4, we show general features of numerical 
solutions taking the case a — oo as an example. The fea- 
tures shown in Figs. 3 and 4 are found also in the cases 
of finite a. In Fig. 3, we show S and 2(1 — a) at ori- 
gin as a function of time t for black hole formation case 
(Cs = 15, dotted lines) and no formation case (Cs = 13, 
solid lines), respectively. The reason why we plot 2(1 — a) 
for a is that it must coincide with 5 for t — * 0. Fig. 3 
clearly shows that for t <C ro, (a) S and a — 1 are pro- 
portional to t, and (b) S agrees with 2(1 — a). Hence, 
the numerical solution reproduces the analytic solution 
for t <C ro accurately. It is found that the difference be- 
tween 8 and 2(1 — a) gradually becomes appreciable at 
t/ro ~ 0.1 and becomes as large as unity when t/ro ~ 1. 
This is reasonable because the density fluctuation ampli- 
tude can become non-linear only after its length scale is 
smaller than the Hubble horizon scale. 

In Fig. 4, we also show ipo as a function of t. As we 
found in Sec. Ill, it does not change so much for t < ro, 
but slightly decreases with time evolution. This is due 
to the effect of the 0(e 2 ) part in ipo- In some previous 
works, ipo at t ~ ro is assumed to be equal to ipo at t — 0. 
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But as shown here, the assumption is strictly speaking 
not correct. 

In Figs. 5, we show the mass of the apparent horizon 
(Mah/^o) as a function of time (i/ r o) m the black hole 
formation cases for a — 2, 3, 5 and oo. We note that in 
all the simulations, computation is terminated when it is 
difficult to keep the numerical accuracy near the apparent 
horizon. Although the simulations had to be ended be- 
fore we could draw a definite conclusion, our results given 
in the figures strongly suggests that Mah approaches an 
asymptotic value without increasing forever. This is con- 
sistent with a previous result |L8 11 1. (We note that by 
comparing it with MKodama, it is found that Mah shown 
here is accurate to within a few percent. ) 

The formation epoch of a black hole is highly depen- 
dent on the initial density profile. In the case a = oo, the 
formation epoch is t > lOOro, but in the other cases, it is 
earlier for smaller a. This is natural because for smaller 
a, the wavenumber of the peak in the spectrum is larger. 
It is found that Mah can be > 8tq in the case a = oo, 
but it is at most ~ l.lr in the case a = 2. These facts 
also indicate that the formation epoch and the maximum 
mass of black holes are highly dependent on the initial 
density profile. We also note that even when we pay at- 
tention only to a particular spectrum shape, Mah can 
vary in a large range depending on Cs (or ipo at t = 0). 
Since we do not pursue the simulation in which a black 
hole in the limit A/ah = is formed, we cannot draw a 
strong conclusion, but it seems possible to have a black 
hole whose mass is much smaller than r near the thresh- 
old of black hole formation. 

In Figs. 6, we plot Mah as a function of ipo at t = 
for the cases a — 2, 5 and oo as examples. We find that 
the threshold of ipo at t — for formation of black hole is 
quite different among all the models. In the case a = oo, 
the threshold value is high {ipo{t — 0) ~ 1.79). On the 
other hand, it is not so large in other cases, and smaller 
for smaller a. In Fig. 7, we plot the threshold line on the 
(a,il)o(t — 0))-plane. The reason of this property is sim- 
ple: In the case in which the density peak is surrounded 
by a high density region, it is forced to be expanded by 
the surrounding region more strongly and is prevented 
from collapsing. As a consequence, the density peak sur- 
rounded by a higher density region is more difficult to 
form a black hole, while the density peak surrounded by 
a lower density region can more easily form a black hole. 
Probably, the density peak surrounded by a zero density 
region can most easily form a black hole. Motivated by 
this observation, we also perform simulations in the flat 
spacetime from the time symmetric initial condition with 
the initial density fluctuation profile, 



2 

fi ( r 

5 ■ ib — C$ exp — 7y 
V ri 



(4.27) 



In this case, the threshold value for ipo at t = is about 
1.428. This value can be approximately regarded as the 
smallest value of ipo at t = for the formation of a black 
hole. 



As another useful criterion, we point out an approx- 
imate measure for determining the formation of black 
holes in the special case of spherical symmetry. First, 
by using the Kodama mass, we define an excess mass at 
each radius at t — as 



SM K (r)= Mr 



M f 



Airp a" 



where Mp denotes the mass in the case of a non- 
perturbed universe and we have used the relations which 
hold at t = such as u = 0(t 2 ), a = 1 + 0(t), and 
d t t/j = 0(t) to neglect the terms higher order in t. Then, 
we define a compaction function at each radius as 



C(r) 



SM K 
rip(r) 2 t 



(4.29) 



Here, we note that this compaction can be defined at 
t = 0, because it approaches a time- in depen dent function 
in the limit t — * 0. In Fig. 8, we show (4.29) as a function 
of r for some of filled circles on the thick dotted line in 
Fig. 7. It is found that the maximum value C(r) max is 
about 0.4 irrespective of a. We note that for each a, 
if the maximum value is large than the value shown in 
Fig. 8, a black hole is always formed, while if not, it is not. 
Thus, the measure presented here will be helpful as an 
approximate criterion to know whether a spherical black 
hole is formed only from the initial data at t = 0, and if 
we use C(r) max as a parameter instead of ipo{t = 0), we 
can approximately neglect the dependence on a. 



V. SUMMARY 

We have presented a formulation for numerical rela- 
tivity in the cosmological background by which we can 
perform a numerical simulation of primordial black hole 
formation with the initial data that can be analytically 
given. Namely, we have formulated the Einstein and hy- 
drodynamic equations in the constant mean curvature 
time-slicing in a way suitable both for obtaining the an- 
alytic solution of a perturbation while its length scale is 
well over the Hubble horizon scale and for performing a 
numerical simulation until the formation of a black hole. 
As a result, it becomes possible to investigate the pri- 
mordial black hole formation from a very early phase of 
the universe just after inflation up to the formation of 
black holes in a continuous manner, without changing 
the gauge condition or numerical matching of the initial 
data. 

By using our formalism, we have carried out a numeri- 
cal simulation of the black hole formation in a spherically 
symmetric, radiation-dominated universe, starting from 
a realistic initial data which are given by the analytical 
solution of a superhorizon scale perturbation. In this pa- 
per, we have considered the initial conditions which are 
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specified by two parameters; one characterizes the ampli- 
tude of the density (or metric) fluctuation and the other 
the shape of the density profile. It is found that the for- 
mation criterion is moderately dependent on the initial 
profile of the density fluctuation. In the case when the 
density peak is surrounded by a flat Friedmann universe, 
the threshold value of ipo at t = for forming a black 
hole is very large ~ 1.8, while when surrounded by a low 
density region, it may be as small as ~ 1.4. This prop- 
erty suggests that the formation of primordial black holes 
may not be determined by a local criterion: Even when 
there is a density fluctuation of a high density contrast, 
it may not efficiently collapse into a black hole if it is in 
a high density region, but it efficiently collapses if it is in 
a low density region. As we have noted, this moderate 
variation of the formation criterion is translated to a very 
large variation in the actual number of primordial black 
holes. Thus, we conclude that the spatial correlation of 
primordial density fluctuations is crucially important in 
studying the formation of primordial black holes. 

In this paper, we have assumed the spherical symmetry 
and restricted our attention to a model which contains 
only two parameters. Even in this case, the formation 
criterion was not so simple. In reality, a primordial black 
hole is formed in a spacetime without any spatial sym- 
metries. In such a case, the anisotropic effect will be im- 
portant in addition to the inhomogeneous effect shown in 
this paper, and it is expected that the formation criterion 
will be much more complicated than the present case. 
Apparently, the next step is to carry out non-spherical 
simulations, and even for such a simulation, the formal- 
ism presented here is perfectly applicable. 
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FIG. 1. The spectrum shapes of the density fluctuation 
(exp(-fc 2 /4) - cxp(-fc 2 CT 2 /4)) are shown for several a. The 
solid line denotes the case for a = oo in which the density 
peak is surrounded by a flat universe, and the dotted lines 
denote the cases for 1.5 < a < 8. 



FIG. 3. S and 2(1 — a) at origin as a function of time t 
for black hole formation case (Cg = 15, dotted lines) and no 
formation case (Cs = 13, solid lines), respectively. 
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FIG. 4. The same as Fig. 3, but for ipo as a function of 
time t. 



FIG. 2. ipo at t = as a function of Cs for a — oo, 2, 3 and 
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FIG. 6. Mah/to as a function of ipo at t = for a — 00, 2 
and 5. 
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FIG. 8. 5M K /rip(r) 2 a as a function of r at t = for the 
critical cases of a — 2, 3, 5 and 00 (i.e., for the initial condition 
denoted by filled circles in Fig. 7). The maximum value is 
~ 0.4 irrespective of a. If the maximum value C(r) max is 
larger than ~ 0.4, a black hole is formed, but if not, it is not 
for all a. 



FIG. 7. Summary of numerical results on formation of 
black holes in a — tpo(t = 0) plane. For initial conditions 
located in the region above the thick dotted line, we find 
that a black hole is formed. Note that the solid line reaches 
tpo(t = 0) ~ 1.79 in the limit a — ► 00. 
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